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NASA Technical Memorandum 78678 

AIRPLANE STABILITY CALCULATIONS WITH A CARD 
PROGRAMMABLE POCKET CALCULATOR 

Windsor L. Sherman 

August 1978 

j 

Please make the following corrections: 

Page 15: Sentence after equation (11) should read as follows: 

Equations (9) and (10) were programmed for the calculator and the- 
program is given in appendix B. 

Page 16: Equation (16) should read as follows: 

b 2 

Re (y) = S + T - — 

3 

Page 24, Last sentence: Change step 49 to step 45. 

Page 25: Step 100 should read as follows: 

STO><9 (ga T /2U ss ) sin 2Yss 

Page 26, Step 105: Change - to + 

Step 141: Change RCL 8 to RCLB 

Page 29: Delete the last sentence. 

Page 49: In column headed "Output," change the values of 83 , a 2 , a-j , 

ao» and a 12 to 

a 3 = 1.3980958 

a 2 = 1.1093007 

= -0.0098076 

a 0 = -0.0211448 

a 12 = 0.0373094 
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Equation (3) should read 


*10 = x e 


'In e 

an 10/ 


1/2 


= x e 


f 1 

an 1 0 1 


A/2 


ISSUED NOVEMBER 1978 


NASA Technical Memorandum 78678 


Airplane Stability Calculations 
With a Card Programmable 
Pocket Calculator 


Windsor L. Sherman 
Langley Research Center 
Hampton, Virginia 


NASA 

National Aeronautics 
and Space Administration 

Scientific and Technical 
Information Office 


1978 



SUMMARY 


Programs are presented for calculating airplane stability characteristics 
with a card programmable pocket calculator. These calculations include eigen- 
values of the characteristic equations of lateral and longitudinal motion as 
well as stability parameters such as the time to damp to one-half amplitude or 
the damping ratio. The effects of wind shear are included. Background infor- 
mation and the equations programmed are given. The programs are written for the 
International System of Units, the dimensional form for the stability deriva- 
tives, and stability axes. In addition to the programs for stability calcula- 
tions, an unusual and short program is included for the Euler transformation of 
coordinates used in airplane motions. The programs have been written for a 
Hewlett Packard HP-67 calculator. However, the use of this calculator does not 
constitute an endorsement of the product by the National Aeronautics and Space 
Administration. 


INTRODUCTION 

Over the past several years, the programmable pocket calculator has devel- 
oped into a highly sophisticated device that has almost computer characteristics. 
Because of its sophistication, the newer models are capable of being programmed 
to make very complicated calculations. Since different logics are used in pro- 
grammable calculators and since the available keyboard instructions vary with 
models of different manufacturers, it is necessary to identify the make and 
model of the calculator for which a program is written. The airplane stability 
programs presented in this paper were written for a Hewlett Packard HP-67 card 
programmable calculator; however, its use and identification in this report does 
not constitute an endorsement of the product by the National Aeronautics and 
Space Administration. 

Programs are given for the calculation of the coefficients of the airplane 
lateral and longitudinal characteristic equations, the eigenvalues, and the 
stability parameters such as the time to damp to one-half amplitude or the 
damping ratio. In addition, a unique coordinate transformation program is 
given for transformations between inertial axes and airplane body axes. This 
program requires very few program steps and may be useful as part of a larger 
program. The equations on which the programs are based are given so that the 
programs can be readily adapted to other calculators that have sufficient pro- 
gram capacity. 

The programs presented herein evolved during the study of wind shear and 
its effect on airplane stability and control. These programs proved useful 
in making stability calculations in this study and should be of use in other 
investigations. 



SYMBOLS 


A 

aspect ratio 


a 0» a l t * 

. .,35 coefficients of characteristic equations 

a l 2» a l 3* 

a-| 4 ,. . . elements 

of longitudinal stability determinant 

b 

wing span 


b 2» b 1 / b 0 

coefficients of resolvent cubic 

C D 

drag coefficient. 

D 

PSui s /2 

C D,o 

drag coefficient 

for Cjj - 0 

c D a 

3C D 

9a 


C L 

lift coefficient. 

L 

PSu| s /2 


lift coefficient 

at zero angle of attack 

c L a 

9 Cl 

9a 


c L<i 

. 9Cl 
" 9a 


C L§ 

3C L 
= 9(5 


Cl 

M X 

rolling-moment coefficient, 

PSbU§ s /2 

C! P 

9Q 

9p 


C ’r 

9Cz 

= 97" 



2 



% 


3 c i 

alT 



3Cj 

33 



3 Ci 
3<t> 


Cm 

Cm,o 

C % 


My 

pitching-moment coefficient, 

P Scul s /2 

total pitching-moment coefficient at zero angle of attack 

3a 


^ ’ is - 

Cin e = a T 


c n 



m 9 


yawing-moment coefficient. 


pSbul s /2 


3fr 

3p 

3C r 

3r 


Cn B a$ 

3C n 

Cn 3 = iT 


3 



_ 9C n 


C n<j) 

* 3<t> 


Gp 

thrust 

coefficient 

Op 

x u 

9Cp 

9u 


C Y 

f y 

side-force coefficient, 

pSUi s /2 


3Cy 


Cy 

*P 

= 3p~ 


Cy 

r r 

3Cy 

= i7~ 


c Yg 

3C Y 

= W 


C Yg 

3C Y 
= 3fT 


c n fC 2 i 

b n » b i 2 

/C30 

» b 1 3 , . . ■ 

^ terms in lateral stability determinant 

c 

mean aerodynamic chord 

D 

drag 


f t 

thrust 


F T,tr 

trim thrust 

f t 

x u 

3 Frp 
3u 


F X' F Y' F 

2 forces 

along X, Y, and Z stability axis 

' X «e 

3 p x 

= ^7 
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3f y 

3*7 

9Fy 

36^ 


g 

I x^ I y/ I z 

T xz 
Im( ) 

k X' k y7l 

k z^ k xzJ 

L 

M X fMyfM Z 


9F Z 

= ^7 

acceleration of gravity 
moments of inertia, stability axes 
product of inertia, stability axes 
imaginary part of complex root 

radii of gyration, stability axes 
lift 

moments about X, Y, and z stability axes' 



9M X 

^7 




M *6 e 

M z 6a 

Mz 6 r 


9m y 

^7 

9m 2 

967 

9m z 

96^ 


n d 

N l/2 


mass 

number of cycles to double amplitude 

number of cycles to damp to one-half amplitude 
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p 

rolling velocity 

R * 

radius 

Re( ) 

real part of complex root 

Re (y) 

real root of resolvent cubic 

r 

yawing velocity 

S 

wing area 

t 

period 

fc D 

time to double amplitude 

t l/2 

time to damp to one-half ampl: 

Uss 

steady-state velocity 

u w 

wind velocity 

u pr 

perturbation velocity 

u w’ 

wind shear gradient 

w w ' 

updraft-downdraft gradient 

X,Y,Z 

stability axes 

x b' Y b* z b 

airplane body axes 

X e / Y e , Z e 

Earth-fixed axes 

*sp'^sp' Z 

s p space axes 

x/y,z 

general variables 

x b'¥b / z b 

body axis coordinates 

x sp/y S p^ z 

S p space axis coordinates 

a pr 

perturbation angle of attack 

a tr 

trim angle of attack 

a l » a 2\ 

a 3 ,a 4 J 

real roots 


3 

Ypr 


6 


sideslip angle 

perturbation flight-path angle 



Y ss 
A 

«a 

«e 

«r 

e l » e 2 


c 

0 tr 

p 

O qi 


a 


u 


steady-state flight-path angle 

logarithmic decrement 

aileron deflection 

elevator deflection 

rudder deflection 

angles 

damping ratio 
trim pitch angle 
atmospheric density 
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u ss w w' 

a w = 

g 

^,0,4) airplane yaw (heading), pitch, and roll angles, respectively 

oo n undamped circular frequency 

Dot over a symbol indicates differentiation with respect to time- 


EQUATIONS PROGRAMMED AND PROGRAM DESCRIPTIONS 

Six programs are presented in this paper- The first three calculate the 
elements of the lateral and longitudinal stability determinants and the coeffi- 
cients of the characteristic equations- In addition, program 3 extracts a real 
root of a fifth-order polynomial when required. Programs 4 and 5 complete the 
root extraction process and calculate the stability parameters- Program 6 
implements the Euler angle transformation by using the polar-rectangular keys 
found on calculators. 

Programs 1, 2, and 3 are written for the International System of Units, 
stability axes (fig. 1), and the dimensional form of the stability derivatives. 
The equations programmed are the linearized form of the equations of motion 
derived in appendix A of reference 1; thus, the effects of wind shear are 
included . 
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In deriving these equations, head winds and updrafts were taken as negative. 
Thus, a positive u w ' will change a head wind into a tail wind, and a positive 

w w ' will change an updraft into a downdraft. The signs of u w ' and w w ' set 

the signs of a u and 0 W ; u w ' is a gradient with altitude and w w ' is a 

gradient along the flight path. 

In writing the programs, the following conventions were used for the labels: 

(1 ) Capital letters (A to E) are program labels 

(2) Lower-case letters (a to e) are subroutine labels 

(3) Numbers (0 to 9) are used for all other labels 

Table I summarizes the programs presented in this paper. The key entries 
given in appendixes A to F are the standard HP-67 key entries given in the 
owner's manual. Check cases for all programs are given in appendix G. 

TABLE I.- SUMMARY OF PROGRAMS 


Program 

Description 

Key entries 
given in 

1 

Calculates the elements of longitudinal stability determi- 
nant and normalized coefficients for characteristic 
equation 

Appendix A 

2 

Calculates the elements of lateral stability determinant 
and starts calculating coefficients of the characteristic 
equation 

Appendix B 

3 

Label A completes calculating coefficients of characteris- 
tic equations of lateral motion; label B calculates a 
real root of a fifth-order polynomial and reduces the 
fifth-order polynomial to a fourth-order one; t-j /2 or 
t D for the real root determined; label B can be used as 
a stand-alone program 

Appendix C 

4 

Uses Ferrari's method to calculate the roots of a fourth- 
order polynomial and can be used as a stand-alone pro- 
gram; will also determine roots of cubic, quadratic, and 
first-order equations 

Appendix D 

5 

Calculates stability parameters such as /2 #■ t D , 

and N-j /2 

Appendix E 

6 

Uses the polar-rectangular transformations of the calcula- 
tor to implement the Euler transformation between space 
and body axes or body and space axes; this method saves 
about 57 program steps when compared with the more usual 
methods of programming 

Appendix F 
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Programs 1 and 2 give solutions from an equilibrium flight condition. 

There are six parameters, U ss , Y S s' a tr' F T tr' °T * and a W ' that must be 
adjusted correctly to obtain the equilibrium flight condition. There are two 
equations to accomplish this adjustment. Programs 1 and 2 were set up in the 
following manner. The parameters U ss , Y S s' 0 T* and a w are specified by the 
user. The program calculates a t r , assuming that F T , tr is 0. For the flight 
condition U ss = 77.12 m/sec, Yss = -0.05236 rad, 0" T = 2.0, and O w = 0.0, the 
error introduced in a tr by this method is 0.00081 rad, which is considered 
acceptable. If it is desired to monitor the calculated value of a tr' insert a 
pause after step 45 of program 1* 

Program 1 

The linearized equation of longitudinal motion is in symbolic form 



The characteristic equation for longitudinal stability is obtained from 
the determinant of the 3x3 matrix and* has the form 



d 

+ ai — + an = 0 

at 


( 2 ) 


where ag to a 4 are given by 

a 4 = a 22 “ a 32 < 3a ) 

a 3 = a l 2 ( a 22 ~ a 32> + ( a 23 “ a 33 “ a 25 a 32 + a 22 a 35) < 3b ) 

a 2 = < a 23 a 35 _ a 25 a 33 " a 26 a 32) + a 12( a 23 “ a 33 “ a 25 a 32 + a 22 a 35> 

+ a l 3 ^ a 31 ~ a 21 ) < 3c ) 

9 



a l = a l 2 ( a 23 a 35 ~ a 25 a 33 “ a 26 a 32) + a 3l( a 13 a 25 “ a 14 a 22> 
- a 2 ] (a-| 3 a 35 - a 14 a 32 ) - a 26 a 33 

a 0 = a 33 ( a 21 a l 4 “ a 26 a 12) + a 31< a 13 a 26 “ a l 4 a 23> 
and a-) 2 , a-] 3 , etc., are given by 

/ v "N 


a l 2 " " 


9 a T / C L ‘ 

sin 2y ss + I C 0fO + —Jk , - 


2U 


ss 


m 


9 •? 

a 13 = c L k l - — < a T sin ^ Yss ~ <?w> 
u ss 

a l 4 = + c m a a tr) k 2 * 

a 21 = c D a k 3 

a 22 ^ c L<i + C L q) k 3 

a 23 = C L a k 3 

a 25 = ~^m0 + ( 'nt ( ^ k 4 

a 26 = -( - : m a k 4 

a 3 -| = g (cos Y ss " °T cos 2 Yss> 

a 32 = " u ss + c L0 k 3 

a 33 = 9 (sin y ss - o T sin 2y ss ) 

a 35 = _c ni0 k 4 


(3d) 

(3e) 


(4) 
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PS^SS PSCUgg PSUgg PSCUgg 

where ki = , ko = , ko = , and k 4 = . In addition to 

m I Y 2m 2 I y 

the foregoing equations, the following equations are needed to calculate the 
values of C L , Cq, and a tr at trim: 


2 mg 

C L = 5 ~< a T sin Y ss - CT w + cos Yss) (5a) 

PSUgg 


C D “ c D,o + 


Cr.2 


TTA 


(5b) 


«tr - 


C L “ C L,0 

C La 


(5c) 


Because large changes in forward speed are encountered in wind shear, the 
effects of the u stability derivatives not normally accounted for are included 
in this program. This was done in the following manner: 


D u 


3d 

3u 



(used in eq. 4) 


3l 

L u = — = C L ki (used in eq. 4) 

3u 


My 


u 


3 My 


3u 



+ Vtr k 2 


(used in eq. 4) 


Equations (3) , (4) , and (5) were programmed to calculate the coefficients 
of the characteristic equation, which is equation (2). The key codes for pro- 
gram 1 are given in appendix A. 

The program destroys the original input data but preserves the coefficients 
of the determinant in the secondary registers. The principal output is the * 
normalized coefficients of the characteristic equation which are stored in Rq, 

R-j , R 2 , and R 3 . 
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Program 2 


The linearized equations of lateral motion with the effects of wind shear 
included are, in symbolic form. 


C 11 — + b 13 
at 


d \ d 

- .b 14 -*b 15 


/d \ 2 d 

b 43 — + b 16 — + b 17 

Vat/ dt 


C 21 — + b 22 
dt 


c 30 77 + b 31 
at 


/d \ d d 

b 42 77 + b 23 77 b 32 77 + b 33 


dt 

/ 

d \ 2 d 

TT + b 24 — 


dt 


dt 


dt 


dt 


b 34 — + b 35 
dt 




r 


* 



* 



6 



_ 



F Y{ a F Y(S r 


M X{ a M x<s r 


M Z6 a M z<5 r 


( 6 ) 


and are based on equations (A5) , (A8) , and (A9) of reference 1 . In lineariz- 
ing these equations, it was assumed that no wind gradient existed in the 
Y e derivative in Earth axes. If the wind gradients are zero (i.e., no wind 
shear) , these equations reduce to the standard form of the linearized equa- 
tions of lateral motion that are given in many standard works, such as refer- 
ence 2. The equations are valid in the interval -0.17453 S y ss s 0.17453. 


The characteristic equation is obtained from the 3x3 matrix on the 
left-hand side of equation (6) and has the form 




d 

+ a 1 — + a 0 
dt 


= 0 


(7) 


for o-r ^ 0. 

When a>p = 0, the sq term in equation (7) becomes 0. Equation (7) now 
has one zero root and four finite roots and is solved as a quartic. Program 2 
tests equation (7) and informs the user if a fourth- or fifth-degree polynomial 
is present. The coefficients ag to as are given by 


a 5 = C 30H “ b 43 b 42) < 9a ) 

a 4 = C 11 ( b 42 b 34 “ b 32) “ C 21 ( b 34 " b 43 b 32) + b 31 ” b 43 b 42> 

+ c 30 ( b 24 + b 14 ~ b 43 b 23 “ b 16 b 42) ( 9b ) 
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a 3 = b 13< b 42 b 34 " b 32) ~ b 22< b 34 “ b 43 b 32> + Cl 1 ( b 42 b 35 + b 23 b 34 “ b 24 b 32 

- b 33 ) - C 2 ] ( b 35 + b l 4 b 34 ~ b 43 b 33 “ b l 6 b 32> + b 31 < b 24 + b l 4 “ b 43 b 23 

~ 6 b 42 > + c 30< b 15 + b 14 b 24 ~ b 1 6 b 23 “ b l 7 b 42> < 9c ) 

a 2 = Cg i (b 23 b 3 5 - b 24 b 33 ) + b i 3 (b 42 b 35 + b 23 b 34 - b 24 b 32 - b 33> " c 2l( b 14 b 35 
+ b 15 b 34 _ b 16 b 33 " b 1 7 b 32> + C 3 o( b 15 b 24 “ b l 7 b 23) “ b 22< b 35 + b 1 4 b 34 

- b 43 b 33 “ b l6 b 32) + b 31 ( b 15 + b 1 4 b 24 “ b 16 b 23 “ b 1 7 b 42> ( 9d ) 

a l = b l 3 ( b 23 b 35 " b 24 b 33) _ b 22 ( b l 4 b 35 + b l 5 b 34 - b 16 b 33 " b 17 b 32) 

- C 2 ] { b-j 5b 35 - b 17 b 33 ) + b 3l( b l5 b 24 “ b l 7 b 23? ( 9e > 

a 0 = _b 22 ( b l 5 b 35 - b l 7 b 33> (9f) 

The C and b terms that are used to generate the coefficients of the charac- 
teristic equation are given by 


b] i - 0*0 


b 1 2 C Yp k 5 


b 1 3 


gO 


W 


n2 

U SS 


Y ss 


J ss 


b 14 = -c Zp k 6 

b 1 5 = U w'(-Cl r k 6j = Cz<t> k 6 
b 1 6 = -c n p k 7 

b l 7 = u w'(^“ < - : n r k 7^ = C 'n(j ) k 7 
b 21 = “ c Y r k 5 


(10a) 

(10b) 

(10c) 

(10d) 
(1 0e) 
(1 Of) 
(10g) 
(1 Oh) 
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g(°u + ° w ) 


b 2 2 - . sin 2 Y SS 

2uls 

(10i) 

b 23 = ~ c l r k 6 

noj) 

b 24 * ~ c n r k 7 

(10k) 

b 30 ■ 

(101) 

b 31 = ~ c Y^ k 5 

(10m) 

b 32 = -C^k 6 

(lOn) 

b 33 = -c 2g k 6 

(lOo) 

b 34 = ~Cng k 7 

(10p) 

b 35 = “C n gk 7 

(lOq) 

*XZ 

b 42 “ “ “ — 
J X 

(10r) 

' Z XZ 

b 43 “ - ~ — 

1 Z 

(10s) 

cn = bn + b i2 

(10t) 

C21 = 1 + b 21 

(10u) 

C 30 = 1 + b 30 

(10v) 


PSU ss PSbU f s PSbUg S 

where k 5 = , k 6 = , and k 7 = . The trim angle of attack 

2m 2I X 21^ 


14 



was calculated from 


a tr 



sin 2 Y ss 


a„ + cos Y ss ) 1 



( 11 ) 


Equations (9) , (10), and (11) were programmed for the calculator and the program 
is given in appendix B. The stability derivatives C n g, and c Yg 


have been included in this program. 


The derivatives 


and 



are always 


calculated when wind shears are included. This program calculates all b and 
C coefficients in the determinant and starts calculating the coefficients of 
the characteristic equation. 


Program 3 

Program 3 completes the calculation of the coefficients of the character- 
istic equation and tests the contents of register 4 to determine if the equa- 
tion is a quartic or a quintic. If it is a quartic, the number 4 is displayed 
and the calculator stops. If it is a quintic, the number 5 is displayed and 
the calculator stops. Label B of this program calculates a real root of the 
quintic equation by using the secant method. The initial guess for the root 
is obtained by dividing the coefficient a 4 by 5; this operation is done in 
the program. 

Even with an estimate of the root, however, two estimated points are 
required to start the secant method. These points are obtained by either add- 
ing or subtracting 0.08 from 34/5. The number 0.08 has proved satisfactory 
for several different fifth-order polynomials. However, the secant method is 
sensitive to this number and changes may be, necessary. Subsequent estimates of 
the root were calculated from 


x i+1 


x i - x i-l 

*i - f(Xi) 

f( Xi ) - f( Xi ^) 


( 12 ) 


where x^ is always the present value and f(x^) is the value of the function 
being used for x = x^. Synthetic division is used to determine when a root had 
been found. The fifth-order polynomial is then reduced to a quartic for pro- 
cessing by program 4. During the iteration process, the calculator pauses to 
display the value of the characteristic polynomial so that convergence can be 
monitored. When the display shows zero, the root has been found. When the 
root has been found, the calculator will stop and display 8. The root is in 
the Y stack register and the time to damp to one-half amplitude or time to 
double amplitude is in the Z register. A negative number in the Z register 
means that the value given is the time to double amplitude. The number of 
iterations required to extract the root is in the T stack register. 



The test used for the determination of a root is that the polynomial must 
be zero to the number of digits in the calculator display; thus, the test for 
a root assures its accuracy. 


Program 4 

A quartic equation is the highest order polynomial for which an explicit 
analytical solution for the root exists. Ferrari's method (refs. 3 and 4) and 
appendix H was used to obtain the roots of the quartic from the characteristic 
equation. The general form of the quartic equation is 

a 4 X^ + a 3 X 3 + a 2 x 2 + a-|X + ag = 0 (13) 

The first step in applying Ferrari's method is to normalize equation (13) so 
that a 4 == 1 . The determination of a real root of the following resolvent 
cubic is the next step: 

y 3 + b 2 y 2 + b-jy + bp = 0 (14) 

The coefficients of equation (14) are given by 


b 2 = -a 2 
b l = a l a 3 “ 4a 0 

b 0 = a o( 4a 2 “ a 3 2 ' " 




a l 




(15) 


and the root Re(y) is obtained by 


b 2 


Re (y) = S + T = 

3 


(f > 0) 


(16) 


or 


Re(y) = 2(R 3 + f)l/3 cos 


-(tan -1 £ 
3 V R 


b 

2 


(f ^ 0) (17) 


where 


Q = (3b-] - b 2 ^) j 9 
R = (9b 2 b] - 27b 0 - 2b 2 3 )/54 
f = R 2 + Q 3 


(18a) 

(18b) 

(18c) 
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S = (r + \/f)V3 


(18d) 


T = (R - \|f)V3 (1 8e) 

The root Re(y) is any root of the resolvent cubic, equation (14); this pro- 
gram is written to calculate the largest real root of equation (14). Once 
Re(y) is known, the roots of the quartic are obtained by solving the following 
two quadratic equations: 


z 2 + (A + C) z + (B + D) = 0 
z 2 + (A - C) z + (3-D) =0 


(19) 


where 

a 3 

A = — 
2 


Re (y) 

B = 

2 


D = 



C = 



( 20 ) 


(D ? 0) 


C 



a 2 + Re(y) 


(D = 0) J 


Equations (15) to (20) and a quadratic solution routine were programmed to 
obtain the roots of a quartic equation. The key codes for program 4 are given 
in appendix D. 

Because f and D are tested to determine program direction, special pro- 
gramming is required both to insure that nonsignificant digits do not influence 
the test and to protect against the small difference of large numbers. The 
expressions for f and D were written as 
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for programming. In each case, the quantity in the parenthesis was rounded to 
the calculator display and then tested. Special routines were added to protect 
against R and B being equal to 0 . The introduction of rounding will intro- 
duce some error if a significant number is truncated. As the rounding is con- 
trolled by the number of decimal digits in the calculator display, there is 
flexibility in the amount of rounding introduced. Experience with a set of 
20 test equations indicates that a display of 7 digits is satisfactory for most 
cases. 

The roots of the quartic are stored in registers R*j , R2, S] , and S2. The 
root indicator (-1.0 for complex roots and 0.0 for real roots) is stored in 
registers Rg and Sg. If the roots are complex, the real part is stored in 
register 1 and the imaginary part in register 2. 

This program is a general program for the roots of a quartic equation and 
may be used as a stand-alone program if the coefficients of the quartic are 
stored in the following locations: 

a3 in register Rg 

a 2 in register R-| 

a-| in register R2 

ag in register R3 

In addition, this program may be used to solve for the roots of lower order 
equations. For the cubic where the equation has the form 

a3X^ + a2X^ + a-jx -1- ag = 0 , 33 = 1 . 0 

the equation is multiplied by x so that it is converted to a quartic with a 
zero root and the coefficients are stored as follows: 

a2 in register Rg 

a-j in register R-] 

ag in register R2 

0.0 in register R3 

Quadratic and first-order equations may be solved in a similar manner by 
multiplying through by x 2 or x^, respectively. 
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Program 5 


Program 5 calculates the stability parameters (ref. 5, p. 61), such as 
the time to damp to one-half amplitude or the damping ratio. The equations 
programmed are given as follows: 

Time to damp to one-half amplitude t*| or time to double amplitude t D : 

0.693 

t V 2 or = (2D 


Per iod : 


t 


2tt 

lm( ) 


r 


Z 'I 

caJ 


(22) 


Number of cycles to damp to one-half amplitude N-j/ 2 
amplitude N D : 


lm( ) 

^1/2 or Nn = -0.110 

7 Re ( ) 


or time to double 


C ) ~ 


(23) 


Logarithmic decrement: 


0.693 

A = 

N] /2 or N d 

Undamped circular frequency: 


w n = t (Re ( ) ) i 2 + (lm( ))2] 1/2 
Damping ratio: 


(24) 


(25) 


C 


Re ( 


Wr 


(26) 


If A, t, or N is negative, unstable conditions are indicated. For instance, 
0.693 

if is negative, the time calculated is for doubling the amplitude. 

Re ( ) 

i The key entries for this program are given in appendix E and the storage 
at the end of this program contains all the calculated information concerning 
airplane stability. For real roots, only the time to damp to one-half amplitude 
or the time to double amplitude is calculated. This program may be used as a 

stand-alone program. 


19 



Program 6 


Program 6 uses the polar -rectangular keys of the calculator to implement 
the Euler transformation used in rigid-body rotation. The transformation pro- 
grammed is the transformation that is frequently used in aeronautics 

(fig. 1). The use of the polar -rectangular keys permits a short program for 
this type of transformation. 


The transformation scheme is illustrated through the use of a two- 
dimensional transformation. The coordinates of a point p(x,y) in the 
xy axis system are given in the x'y' axis system, which is rotated through 
the angle t-\ with respect to the xy axis system by 


x' = x cos Ei + y sin £■) 
y' = -x sin e-j + y cos £-| 


(27) 


The polar coordinates of p(x,y) are R^, £2 in the xy axis system, where 

R* = (x 2 + y 2 )^^ and £ 2 = tan - ^ — , and are R^, (£2 - e i) in the x'y' axis 

x 

system. The x'y' axis system coordinates are now given by 


x ' = R^ cos (£2 - e l ) 


y' = R* sin (£2 - £() 


(28) 


If equation (28) is expanded (x is substituted for R* cos £ 2 and Y is 
substituted for R* sin £ 2 ) , equation (27) results and- shows that the same 
transformation is taking place. This result leads to a program for a two- 
dimensional transformation. It is assumed that y is stored in the Y stack 
register, x is stored in the X stack register, and £-| is stored in regis- 
ter Rn. The program is as follows: 


+P 


x+y 
RCL n 


x+y 

-*-R 


This program gives x' and y' in 6 steps instead of the usual 18 steps. This 
two-dimensional program is completely general. If this two-dimensional trans- 
formation program is used in conjunction with a bookkeeping program, three- 
dimensional transformations may be made. In reference 6 (pp. 272 to 275) , a 
method is given that simplifies the bookkeeping problem. A program for one of 
-the three-dimensional Euler transformations used in aeronautics is given in 


20 



appendix F. This program is for transformations between two right-hand axes 
systems (fig. 1) in which the Z-axis is positive downwards. The first rota- 
tion is through the angle ip about the Z S p-axis; the second is through the 
angle 0 about the Y^p-axis; and the third is through the angle cj) about the 
X^-axis. The angles y, 0, and <j) are the airplane heading, pitch, and roll 
angles, respectively. The program presented in appendix F is a specialized pro 
gram because, in three-dimensional transformations, the order in which the rota 
tion angles are taken and the axes about which the rotations take place vary 
from one transformation to another. Similar programs may be written for other 
three-dimensional transformations by changing the bookkeeping part of program 6 
Subroutines B and C would not be changed. 

The advantages of using the polar -rectangular keys in program 6 for three- 
dimensional transformations are not apparent unless program 6 is compared with 
a program that uses the traditional approach of calculating the direction 
cosines and then using them to make the transformation. By using direction 
cosines, a reasonably efficient program for the \p,0,<}) transformation dis- 
cussed in this section takes 124 program steps and 20 storage registers, com- 
pared with 67 steps and 10 storage registers for the polar -rectangular method 
of this paper. The impact is even more apparent if both the polar -rectangular 
(P+R) and the direction-cosine (D-C) methods are considered as subprograms to 
a main program. Take the following example: 

A vector has been computed and its components are stored in three consec- 
utive registers. The angles ty, 0, and <p have also been calculated and are 
stored in consecutive registers. It is desired to transform the calculated 
vector components to a new coordinate system rotated from the original by the 
angles 6, and <{>. 

Table II summarizes the manner in which the two programs would merge with 
the main program. Storage for the original vector components and the angles is 
not counted. 

TABLE IX.- COMPARISON OF THREE-DIMENSIONAL TRANSFORMATIONS WHEN USED IN PROGRAM 


Programming considerations 

Space 

to body 

Body to space 

Two 

way 

P*R method 

D-C method 

P+R method 

D-C method 

P+R method 

D-C method 

Program steps in transformation 

26 

83 

27 

92 

53 

118 

Registers used 

— 

13 

— 

13 

— 

13 

I register 

Used 

Used 

Used 

Used 

Used 

Used 

Storage for new computation 

3 

3 

3 

3 

6 

6 

Total program steps used 3 

26 

83 

27 

92 

53 

118 

Total registers used 3 

3 

16 

3 

16 

6 

19 

Steps available for main program 

198 

141 

197 

132 

1 71 

106 

Registers available for main 







program 

22 

9 

22 

9 

19 

6 


a I register is not counted. 21 



Analysis of the data presented in table II shows the economics of using 
the P*R method in programs. In addition, a calculator with only 49 program 
steps can be programmed one way using the P*R method, and a calculator with 
98 steps can handle the two-way P*R transformation. For the D-C method, the 
smallest programmable calculator that can handle a one-way transformation is 
one with 98 steps. The two-way transformation will not fit on a 98-step 
calculator . 


USE OF PROGRAMS 1 TO 6 

After a program has been keyed into the calculator, the program switch 
should be set to run. Set the display and trig modes, switch back to program, 
and record program.^ The display and trig mode status are now recorded on the 
magnetic card and the calculator will be set to the indicated status conditions 
whenever the program is read in. The display and trig mode status are given for 
each program in the* appendixes. 

Appendix G contains the check case for the programs in appendixes A to F. 

To make longitudinal stability calculations, use the following procedure: 

(1 ) Enter program 1 (appendix A) 

(2) Enter data as shown on storage map 
Push A 

At stop, coefficients of characteristics equation have been calculated 

(3) Enter program 4 (appendix D) 

Push A 

At stop, roots of characteristic have been determined 

(4) Enter program 5 (appendix E) 

Push A 

At stop, complete set of longitudinal stability data is stored as 
indicated on storage map 

To make lateral stability calculations, use the following procedure: 

(1 ) Enter program 2 (appendix B) 

(2) Enter data as shown on storage map 
Push A 

(3) At stop, enter program 3 (appendix C) 

Push A 

If 4 is displayed at stop, go to step 4 
If 5 is displayed, push B 

When 8 is displayed, the real root, the time to damp to one-half ampli- 
tude or the time to double amplitude, and the number of iterations are 
stored in the stacks 
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(4) Enter program 4 (appendix D) 

Push A 

At stop, roots of quartic have been calculated 

(5) Enter program 5 (appendix E) 

Push A 

At stop, a complete set of lateral stability data has been calculated. 
Data relating to quartic is stored in calculator 

Programs 4 and 5 may be used as stand-alone programs. 

Program 6 may be used in several different ways. To transform from space 
axes (X s p,Y S p, Z s p) to airplane axes (X 5 , Z^) , use the following procedure: 

(1) Enter z sp'Ysp' x sp stac ^ i n or< 3er given 

Push A 

(2) Enter in stack in order given 

Push B 

Push C to make the transformation 

At stop, airplane axis coordinates x b'Yb' z b are stored in registers 
R$, R 7 , and Rg/ respectively 

To transform from airplane axes (X^, Y^, Z^) to space axes (X S p, Y S p, Z S p) , use 
the following procedure: 

(1) Enter z b'Yb' x b stack in order given 

Push A 

(2) Enter in stack in order given 

Push B 

Push D to make the transformation 

At stop, space coordinates x sp'Ysp' z sp are stored in registers R 5 , 

R 7 , and Rq , respectively 


Langley Research Center 

National Aeronautics and Space Administration 
Hampton, VA 23665 
April 28, 1978 
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APPENDIX A 


PROGRAM 1 - LONGITUDINAL AIRPLANE STABILITY 

Program 1 uses the basic physical and aerodynamic data of an airplane 
to calculate the coefficients of the characteristic equation of longitudinal 
motion. This program calculates normalized coefficients for the characteristic 
equation. For running this program, the calculator should be set to calculate 
in radians and to display seven decimal digits. A pause inserted after step v 4! 
causes the calculated value of a tr to be displayed. 
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001 

LBLA 

Calculate elements 


RCL9 

of determinant 


STO+3 

RCLO 

X 2 

O rp 


CHS 

STOr (i) 

-c/ky 2 


RCL 8 

RCL7 


010 

SIN 

X 2 



X 

RCL9 



STOO 

RCL7 

cos 

ST09 

+ 

sin 2 Y ss " 

020 

RCL5 



X 

g (a T sin 2 Yss " 0 w 


RCLB 

RCL2 

X 

RCL 6 

X 

RCL1 

+ cos Yss) 


S TOE 

pSU ss /m 

030 

RCL 6 

X 



2 ' 

2 


STOx3 
ST OB 

P 


RCL4 

P*S 

C L 


ST06 

stored in R 6 

040 

X->Y 


ST05 

RCLD 

RCLl 

0 ^ stored in S 5 


f 

- o t r 


X 

STO+9 

RCL5 

X 2 


050 

RCLC 



7T 



STO+4 

RCL2 

c D ,o + c l 2 Aa 


STO+3 

RCL 8 

C L 0 + Clj a 


ST 0+7 
RCLE 

c % + ^ 

060 

STOx 4 
STO*5 
RCLA 

X 



STOx 9 
RCLB 
STOx (i) 

a l 4 


STOxO 

a 21 


STOxl 
STOx 2 

a 23 

07 0 

STO* 3 
RCLA 

a 22 


STO ><6 

a 26 


STOx7 

a 25 


STOx 8 

P*S 

RCL9 

RCL7 

2 

X 

a 35 

08 0 

COS 

RCL 8 

STOx9 

X 

STOA 

RCL7 

SIN 

STOx9 

RCL9 


090 

2 

X 

RCL5 



STOx (i) 

X 

a 31 


STOB 

RCL5 

RCL 6 

a 33 

100 

STOx 9 

9°T sin 2 Y ss / 


25 



APPENDIX A 


no 


120 


130 


140 


150 


RCLO 

x 

RCL9 

RCL3 


g(a T sin 2 Y S s - ^)/U ss 


RCLA 

RCLO 

RCL5 


RCL 6 


P*S 


STO- 2 

a 32 

R^ 


STO- 4 

a l 2 

R+ 


STO- 5 

a 1 3 


1 

0 


STOI 

RCL3 

RCL2 

STOC 

RCL4 

x 

RCLl 

RCLB 


+ 

STO(i) 
160 ISZ 

RCLE 
RCL4 

x 

RCL 6 

RCLB 

x 

RCL5 

RCL7 

170 x 

RCL9 

RCL3 

x 

RCLA 


RCL2 + 

RCL7 RCL5 

x RCL 8 

180 x 

RCL3 RCL9 

RCL 8 ' RCL2 

x x 

+ 

STOD RCLO 

+ x 

STO ( i) a 3 - a I 

ISZ STO ( i) 

RCLl RCLO 

RCL 8 190 RCL9 

x x 

RCL7 RCL 6 

RCL 8 RCL4 

x x 


RCL 6 

RCL2 


STOE 

RCLD 

RCL4 

x 

+ 


RCLB 

x 

RCL5 
RCL 6 
200 x 

RCL9 

RCLl 

x 
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RCLA 

x 

+ ao 

P*S 
ST03 
21 0 RCLC 
STOf 0 
STOf 1 
STOf 2 
STOf 3 
RTN 
' R/S 
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Storage Map for Program 1 


Address 

Register 

Input storage 

Output storage 

0 

R o 

ky 

3a 3 

1 

R i 

m 

a 2 

2 

r 2 

P 

a l 

3 

r 3 

Op 

A u 

a 0 

4 

r 4 

C m a 


5 

r 5 

g 


6 

r 6 

u ss 


7 

r 7 

Yss 


8 

r 8 

°u 


9 

r 9 



10 

s 0 

° D a 

a 21 

11 

Si 


a 23 

12 

s 2 

C L0 

a 32 

13 

s 3 

° L a 

a 22 

14 

s 4 

c D f O 

a l 2 

15 

s 5 

0.0 

a l 3 

16 

s 6 

0.0 

a 26 

17 

s 7 


a 25 

18 

s 8 

c mQ 

a 35 

19 

s 9 

c m,o 

a l 4 

20 

r A 

c 

a 31 

21 

r b 

S 

a 33 

22 


A 


23 

r d 

c L,o 


24 

r E 

0.0 


25 

I 

20 



a These are the normalized coefficients of the quartic; 

84 = 1.00. 


thus , 



APPENDIX B 


PROGRAM 2 - LATERAL AIRPLANE STABILITY 

Program 2 uses the basic physical and aerodynamic data of an airplane to 
generate the coefficients of the lateral stability determinant. After complet- 
ing the calculation of these coefficients, the program starts but does not fin- 
„ ish calculating the coefficients of the characteristic equation for lateral 
motion. For running this program, the calculator should be set to calculate in 
radians and to display seven decimal digits. A pause inserted after step 44 
causes the calculated value of a tr to be displayed. 
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001 


010 


020 


030 


040 


050 


LB LA 



RCLO 


RCLl 



X 

<3°w 

RCL2 



- 


RCL3 



RCL2 

Us s 

X 



RCLO 


X 



RCLl 


2 



T 


T 



RCL6 


RCLA 



STO x O 


T 

pSU ss /2m 

060 

X 


CHS 



■f 


ST03 



2 


RCLl 



V 


STOtO 

g/U ss in Ro 


RCLC 


X 



2 


STOx4 

2 

pSbU ss /2m 


X 


RCL4 



sin 


STOE 



X 


RCL5 



STOA 


X 2 


07 0 

R+ 


ST05 

2 


STOC 


ST Of 4 

P SbU ss /2l x 


RC L3 


RCL9 



GSBa 

b l 2 i 

X 2 



DSZ 


ST 09 

2 


GSBa 

b 21 

STOt ( i ) 

pSbU ss /2l z 


1 


RCL8 



STO+ { i ) 


X<0 



X+Y 


SF2 



DSZ 


X 2 


08 0 

GSBa 

b 30 

F?2 



X-^Y 


CHS 



STO+ ( i ) 

c 30 

CHS 



X+Y 


ENT 



GSBa 


ENT 



RCL4 


RCL9 



GSBa 

b 32 

T 

I xz/ I z 


GSBa 

b 23 

ST 09 



RCL(i) 


x+y 



ST06 


RCL5 


09 0 

Rl 


f ’ 



GSBa 

b 1 4 

ST 08 

I xz/ I x 


GSBa 

b 33 

RCLO 



RCLE 


RCL7 



GSBa 

b 34 

RCLl 



GSBa 

b 1 6 

T 



GSBa 

b 35 

X 



GSBa 

b 24 

ST02 



RCL(i) 


RCLC 



ST 07 


COS 






9 


U 


ss 


cos 'Y ss 


C21 
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100 RCLO 
STOx6 

STO*7 

2 

3 

STOI 

RCL8 


P>S 

RCL3 

x 

110 RCL7 

STOE 

RCL6 

RCLl 

x 

RCLO 

RCL4 

x 


1 20 RCL6 
RCL3 
x 

RCLO 

RCL7 

x 


RCL4 

RCLl 
130 P*S 

RCL8 
x 
+ 

GSBb 

ST02 

Rl 

STOI 

Rl 

GSBb 
140 ST03 

Rl 

ST 0+2 

RCLE 

GSBb 

STO+1 

R* 


Calculate coeffi- 

150 

STOO 

RCL9 

P>S 

RCL7 

X 

CHS 

RCL3 

Secondary called 

cient of charac- 
teristic equation 
Secondary called 


+ 

RCLl 
RCL5 
RCL 3 

X 

+ 

b 34 - a 32 a 34 

b 42 b 34 ~ b 32 

160 

RCL2 

RCL7 

X 



RCL4 




P*S 

RCL9 

y 

Primary called 

b 23 b 35 - b 24 b 33 



b 35 + b 34 b l 4 



DSZ 

“ b l 6 b 32 ~ b 33 b 43 


170 

DSZ 

GSBb 

Complete calcula- 



ST 0-2 

tions and store 



Rf 

ST O-l 
Rl 

GSBb 

STO-1 

R4- 

terms 

Primary called 

180 

STO-O 

RTN 

End of program 



LBLa 

Subroutines for 

b 23 b 34 " b 24 b 32 


DSZ 

calculation of 

- b 33 + b42b" 35 


STO* (i) 

coefficients 

Complete calcula- 


RTN 


tions and store 

190 

LBLb 

ENTt 

ENTf 

RCL(i) 

X 

X+Y 

DSZ 

RCL(i) 

X 

ISZ 

RTN 

R/S 
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Storage Map for Program 2 


Address 

Register 

Initial storage 

End of program 

0 

R 0 

g 

(a) 

1 

R 1 

Uss 

(a) 

2 

R 2 

P 

(a) 

3 

r 3. 

S 

(a) 

4 

r 4 

b 

(a) 

5 

r 5 

k X 


6 

r 6 

a u 

b 15 

7 

r 7 

. ff W 

b 1 7 

8 

r 8 

bk xz 

b 42 

9 

r 9 

k z 

b 43 

10 

s 0 

Cfl r 

b 24 

n 

Si 


b 35 

12 

S 2 

° n p 

b 1 6 

13 

S 3 

C ng 

b 34 

14 

s 4 

C lQ 

b 33 

15 

s 5 

Ci 

Z P 

b 1 4 

16 

s 6 

C 2r 

b 23 

17 

s 7 

c 23 

b 32 

18 

s 8 

C Y3 

b 31 

19 

S 9 

c Y3 

c 30 

20 

r A 

m 

b 22 

21 

r b 

C Y r 

C 21 

22 


Yss 

b 1 3 

23 

r d 

Cy 

Y P 

C 11 

24 

r e 

0.0 


25 

I 

24 



of 


a Registers Rg to R4 contain the partially calculated coefficients 
the characteristic equation* 

^If ^xz imaginary, enter as a negative number. 
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PROGRAM 3 - LATERAL AIRPLANE STABILITY (Concluded) 

Program 3 completes the calculation of the coefficients of the character- 
istic equation of lateral motion that was started in program 2. The program 
then determines if the characteristic equation is a quartic or a quintic. If 
it is a quartic, a 4 is displayed and the program stops. Program 4 is then 
used to obtain the roots of the quartic. If the characteristic equation is a 
quintic, a 5 is displayed and the program continues on to extract the real root 
of the quintic and then calculates the time to damp to one-half amplitude or 
the time to double amplitude. This program uses the storage that existed at 
the end of program 2. This program calculates normalized coefficients for the 
quartic and the quintic. For running this program, the calculator should be 
set to calculate in radians and to display seven decimal digits. 
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001 LB LA 

RCL7 
RCL 6 
P+S 
RCL3 
x 

X+Y 

RCL7 

x 

01 0 

RCL2 

RCL4 

x 

RCL5 

RCLl 

x 

+ 

RCL4 
020 RCLl 
P+S 
RCL 6 


RCLO 

P+S 

RCL 6 

Secondary called x 

X->-Y 

RCL7 

x 

060 

DSZ 

DSZ 

GSBa 

STO+3 

R+ 

STO+2 

Rl 

GSBa 

STO+1 

b 1 5 b 34 “ b 32 b 1 7 070 

“ b 16 b 33 + b 1 4 b 35 ST0 ° 

1 

Primary called RCL 8 

RCL9 


x 

X+Y 

RCL7 

' x 

GSBa 

CHS 

030 ST04 

R4- 

STO-3 

R-t- 

GSBa 

STO-3 

Rl 

STO-2 

RCL 8 

RCL9 

040 P+S 

RCL 6 

x 

CHS 

X+Y 

RCL2 

x 

RCL5 

+ 

050 RCLO 

+ 

RCL 6 


x 

RCL 6 
RCL7 
RCL 8 

x 

P+S 
RCL2 
RCL 6 

x 

RCL5 
' RCLO 

x 

090 + 

Secondary called P+S 

GSBa 

STO+2 

Ri 

STO+1 

Rl 

GSBa 

STO+0 

R+ 

bi 4 + b24 — b"j gb42 1 00 STOE 
— b43b23 RCL4 

X/0 

GOTO! 


b 35 b 1 5 " b 33 b 1 7 
Complete and store 080 

calculations 


Primary called 


b 24 b 1 5 - b 23 b 1 7 


1 - b4 2 b 43 


Secondary called 


b 1 5 “ b 1 7 b 42 

- b 1-6 b 23 + b 1 4 b 24 
Primary called 

Complete and store 
calculations 


Determines if equa- 
tion is a quartic 
or a quintic 
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no 


120 


130 


140 


150 


GSBd 



RCL8 


4 

Indicates quartic 


RCL9 


RTN 

Stop for quartic 


ST08 


LBL1 



X*Y 


GSBd 



- 


5 

Indicates quintic 

160 

f 


RTN 

Stop for quintic 


RCL8 


LBLB 

Calculate real root 


X 


RCLO 

of quintic 


- 


STOA 



ST06 


RCLl 

This section 


GOTOO 


STOB 

positions data 


LBLl 

Output routine 

RCL2' 



RCL7 


STOC 



« 


RCL3 



6 


STOD 


170 

9 


RCL4 



3 


STOE 



CHS 


0 

Initialization for 


RCL6 


ST07 

secant method 


T 


FIX 



RCL6 


RCLA 



8 


5 



RTN 

End of program 

• 



LBLa 

Subroutine for 

4 



ENTt 

calculating 

0 


180 

ENTf 

coefficient of 

8 



RCL(i) 

characteristic 

- 



X 

equation 

ST05 



X^Y 


ST06 



DSZ 


GSBb 



RCL(i) 


ST08 



X 


4 



ISZ 


i 



RTN 


6 



LBLb 

Polynomial evalu- 

STO+6 


190 

2 

ation subroutine 

LBLO 

Evaluates polynomial 


0 


1 

and tests for 


STOI 


STO+7 

solution 


1 


GSBb 



GSBc 


ST09 



STOO 


RND 



GSBc 


Pause 

Displays value of 


STOI 


X=0 

polynomial 


GSBc 


GOTOl 

polynomial 


ST02 


RCL6 

Calculates and stores 

200 

GSBc 


RCL5 

new value of X 


ST03 


RCL6 



GSBc 


ST05 

' 


ST04 


X+Y 



RTN 


- 



LBLc 

Synthetic division 





subroutine 
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RCL6 

x 

RCL(i) 

+ 

210 ISZ 

RTN 

LBLd Normalization 

RCLE subroutine 

STOvO 

STOrl 

STOt2 

STOt3 

STOt4 

RTN 

220 R/S 
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Storage Map for Program 3 

(i)Address Register Initial storage 9 End of program * 5 


0 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 
'25 


R 0 

R 1 

R 2 

r 3 

r 4 

r 5 

r 6 

R? 

r 8 

r 9 

s 0 

Si 

s 2 

53 

5 4 

s 5 

s 6 

s 7 

s 8 

S 9 

R A 

r b 

Rc 

r d 

r e 

I 


a 3 

9 2 

a l 

a 0 


b 15 

b 17 

b 42 

b 43 
b 24 
b 35 
b l 6 
b 34 
b 33 
b l 4 
b 23 
b 32 
b 31 
c 30 
b 22 
C 2 1 

b l 3 
Cll 

m use 



a The initial storage is the same as that at end of program 2. The 
partially calculated coefficients of the characteristic equation are 
stored in Rq to R 4 . 

b The end storage is the same for display signals 4 and 8 ; the nor- 
malized coefficients of the quartic are in registers Rq to R 3 . The 
real root of the quintic is in the Y register and the time to damp to 
one-half amplitude or the time to double amplitude is in the Z register 
when 8 is displayed. Pressing R+ moves the real root of the quintic to 
the X register; pressing Rl again moves the time to damp to one-half 
amplitude or the time to double amplitude to the X register. The number 
of iterations required to obtain the root is in the stack T register and 
may be obtained by pressing R I. 
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PROGRAM 4 - ROOTS OF A QUARTIC EQUATION 

Program 4 applies Ferrari's method for the roots of a quartic equation to 
the output of either program 1 or program 3 to determine the remaining eigen- 
values of the characteristic equation of longitudinal or lateral motion. Nor- 
malized coefficients must be used for this program. For running this program, 
the calculator should be set to calculate in radians and to display seven 
decimal digits. 
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001 


010 


020 


030 


040 


050 


LBLA 

Calculate coeffi- 


GOTO! 

4 

cients of resolvent 


LBL0 

ST06 

cubic 


■f 

RCLl 



1 

STOC 



GSBa 

STOx 6 



RCLA 

CHS 



X 

ST04 

b 2 in R 4 

060 

LBLl 

RCLO 



ABS 

STOB 



v 

ST05 



F?2 

X 2 



GOTO 2 

STO -6 



RCL7 

RCL2 



->-P 

STOD 



GSBb 

STO x 5 



2 

X 2 



X 

RCL3 


070 

X-i-Y 

STOE 



3 

STO *6 



T 

4 



COS 

X 



X 

STO-5 

b-) in R 5 


GOTO 3 

x+y 

) 


LBL2 

STO -6 

b 0 in R 6 


RCL7 

3 

' Calculate Q, R, 


X->Y 

STOt 4 

Q 3 , R 2 , and f 


STO-7 

STOr 5 


08 0 

+ 

RCL5 



GSBb 

RCL4 



RCL7 

X 2 

- 


GSBb 

- 

Q 


+ 

X->Y 



LBL3 

STO x 5 



RCL4 

Y x 

Q 3 


- 

RCL4 



RND 

RCL5 



Pause 

X 


090 

ST08 

RCL 6 



ST09 

- 



2 

2 



STOfO 

T 



STOr 2 

RCL4 



STOr 8 

3 



RCLO 

Y x 



ST06 

- 

R 


RCL 8 

ST07 



ST07 

X 2 

R 2 

100 

X 

STOA 



RCL2 

x*o 

• 


- 

GOTOO 

- 


1 

GSBa 



RCL3 


f 

Calculate largest 
real root of 
resolvent cubic 


Re(y) 

Display Re (y) 
Calculate A, B, 
A ± C, B ± D, 
C, D 
A 

B 


C, D fi 0 
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CHS 
RCL8 
X 2 
X/0 
GOTO 4 

no X*Y 
GSBa 
GOTO 5 
LBL4 

GSBa 
RCL8 
. X 2 

x 

LBL5 

120 „ 

STO+7 

STO-8 

F?2 

GOTO 6 

RCLO 

X 2 

RCLl 

RCL9 

130 + 

GOTO 7 
LBL6 

LBL7 
STO+O 
STO-6 
RCL7 
RCLO 
140 GSBc 

RCL8 
RCL6 
P*S 
GSBc 
P+S 
RTN 
LBLa 
+ 

RND 

150 Pause 

X>0 
SF2 
RTN 
LBLb 
X<0 
SF2 


ABS 

3 

1/X 

160 Y x 
F?2 
CHS 
RTN 
LBLc 
X/0 
GOTO 8 
X*Y 

x=o 

GOTO 9 

170 X*Y 

LBL8 

D 2 

B + D t 

B - D CHS 

ST04 

X 2 

X+Y 

ST05 

180 X<0 

GOTOO 


C, D = 0 RCL4 

X<0 

SF2 

X^Y 

F?2 

A + C CHS 

A - C + 

Solve for roots 190 STOt5 

of quartic RCL5 

X+Y 

0 


GOTOl 

LBL0 

ABS 

V 

RCL4 


Subroutine used in 1 

calculation of f 200 CHS 
and D GOTOl 

Displays quantity LBL9 

tested ENTt 

ENT+ 

^ LBLl 

Subroutine for cube STOO 

root of positive or -X- 

negative number R+ 


Guadratic solution 
subroutine 
Protects against 
case 

A±C = B±D = 0 


Calculates -(b/2) 
and (b/2) 2 - C 

Determines if roots 
are real or 
complex 


Solves for real 
roots 


Solves for complex 
roots 


Enters zero roots 


Stores and displays 
roots 
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ST01 
10 -X- 

R4- 

ST02 

-X- 

RTN 



APPENDIX D 


Storage Map for Program 4 


(i) Address 

Register 

Initial storage 3 

End of program* 3 

0 

R o 

a 3 

Root- type indicator 

1 

R i 

a 2 

Re ( ) or a i 

2 

, r 2 

a l 

Im( ) or a 2 

3 

r 3 

a 0 


4 

r 4 



5 

r 5 



6 

r 6 



7 

r 7 



8 

r 8 



9 

r 9 



10 

s 0 


Root- type indicator 

n 

S 1 


Rei ( ) or 03 

12 

s 2 


Imi ( ) or ot 4 

13 

S 3 



14 

s 4 



15 

s 5 



16 

s 6 



17 

s 7 



18 

S8 



19 

s 9 



20 

r A 



21 

r B 



22 

Rc 


a 2 l Coefficients 

23 

r d 


ai / of quartic 

24 

r e 


a oJ 

25 

I 




a Initial storage is provided by output of program 1 or program 3 . 
^The root-type indicator is 0 for real roots and -1 for complex 
roots. The real part of the complex root is stored in R-] or S] and the 
imaginary part in R2 or S2. 
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PROGRAM 5 - STABILITY PARAMETERS 

Program 5 utilizes the eigenvalues computed by program 4 to calculate sta- 
bility parameters, such as the time to damp to one-half amplitude or the damping 
ratio. For running this program, the calculator should be set to calculate in 
radians and to display seven decimal digits. 
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001 LBLA 
RCLO 
RCL1 
RCL2 
CLRREG 
ST02 
R+ 

ST01 

R+ 

010 STOO 

P+S 
RCLO 
RCLl 
RCL2 
CLRREG 
ST02 
R4- 
STOl 
Rl 

020 STOO 

P^S 


6 

9 

3 

CHS 

STOA 

• 

1 

030 1 

1 CHS 

1 ST OB 

IT 
2 

STOI 

x 

STOC 

LBLB 

RCLO 

040 X?0 

SF2 
RCLl 
ABS 
RCL2 
ABS 
+ 

x=o 

GOTOl 

F?2 

050 GOTO 2 

RCLA 
ST03 


Clears registers ST04 

and protects RCLl 

roots and root STOt3 

indicators RCL2 

STO-j-4 


GOTOl 

LBL2 

060 RCLl 
x=o 
GOTO 3 
RCLA 
ST03 
ST04 
X+Y 


STOf3 

ST08 

RCL2 

07 0 X+Y 

Stores constants RCLB 

and initializes I x 

register ST05 


CHS 
STOt4 
LBL3 
RCLC 
ST06 
080 RCL2 
STOt 6 
RCLl 
+ P 
ST07 
CHS 
STOf8 
LBL1 

Determines if roots P*S 

are real or DSZ 

complex 090 GOTOB 

RTN 

Protects against R/S 

zero roots 


Switch for complex 
roots 

Calculates t -| /2 or 
t D for real roots 


Calculate stability 
parameters for 
complex roots' 
Protects against 
zero real part of 
complex root 


ti /2 or t D 


N-j /2 or n D 


A 


t 


to n 


C 

Switch for second 
set of roots and 
program stop 
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Storage Map for Program 5 


Address 

Register 

Initial 

storage 3 

End of program 

0 

R o 

Root- type 

indicator 

Root- type indicator 

1 

R i 

Re( ) 

or a-| 

Re ( ) or 

a l 

2 

*2 

lm( ) 

or a 2 

Im ( ) or 

a 2 

3 

*3 



fc 1/2 or 

t D 

4 

r 4 



A 


5 

r 5 



Ni /2 or 

n d 

6 

r 6 



t 


7 

r 7 



to n 


8 

r 8 





9 

r 9 





10 

so 

Root- type 

indicator 

Root-type indicator 

11 

Si 

Re! ( ) 

or a 2 

Rei ( ) or 

a 3 

1 2 

b 2 

Im! ( ) 

or U 4 

I mi ( ) or 

a 4 

13 

S 3 



tl /2 or 

fc D 

14 

s 4 



A 


15 

s 5 



N 1/2 or 

n d 

16 

S6 



t 


17 

s 7 



% 


18 

s 8 





19 

s 9 





20 

r A 





21 

r b 





22 






23 

r D 





24 

r e 





25 

I 






a The initial storage is the same as the storage at the end of 
program 4. 


For real roots, only the time to damp to one-half amplitude or the time to 
double amplitude is calculated. This quantity is stored in register 3 for the 
root in register 1 and in register 4 for the root in register 2. 
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PROGRAM 6 - EULER TRANSFORMATION FOR AERONAUTICS 

Program 6 is for the standard Euler transformation that is used in aero- 
nautics between inertial axes and airplane axes. The trigonometric mode and 
the number of decimal digits in the display are assigned by the user. 
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001 LB LA 

STOO 

r\ .L 
I\ T 

STOl 

R+ 

ST02 
RTN 
LB LB 
ST03 

01 0 R^ 

ST04 

R^ 

ST05 

RTN 

LBLC 

3 

STOI 
RCL1 
RCLO 
020 GSBb 

RCL2 
GSBb 
X+Y 
ST06 
R+ 

X+Y 

GSBb 

ST07 

R+ 

030 ST08 

RTN 
LB Lb 
+P 
X+Y 
RCL(i) 

X+Y 

“►R 

ISZ 

040 RTN 

LBLD 
5 

STOI 

RCL2 

RCLl 

GSBc 

X+Y 

RCLO 

X+Y 

050 GSBc 

ST08 
R* 


Stores X sp ,Y sp ,Z sp GSBc 

or X b ,Y b ,Z b ST06 

ST07 

RTN 

LBLc 

*>P 

Stores ip, 6, <|> 060 X+Y 

RCL(i) 

CHS 


Transformation 
subroutine X b 
Y b' z b to ^ S P 
Y sp' Z sp 


X+Y 

•>R 

DSZ 

Transforms X S p, RTN 

Y sp t z sp to R/S 

x b 9 Y b 9 2 b 


Transformation sub- 
routine X sp ,Y sp , 
z sp to x b' Y b' z b 


Transforms X b ,Y b ,Z b 
to ^sp^ Y sp^ z sp 



APPENDIX G 


CHECK CASES FOR PROGRAMS 1 TO 6 

This appendix gives check cases for each program given in appendixes A 
to F. Each check case is complete in itself and does not depend on the out- 
put of a previous program. For program 3, two check cases are given - one for 
label A and one for label B. There is no check case given for programs 1 , 2, 
and 3 for a u = a w = 0.0. All check cases are independent of previous results. 
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Check Case for Program 1 


Register 

Input storage 


Output 

R o 

k y = 10.463784 

a 3 

= 1 .3924836 

R 1 

m = 90909.1 

a 2 

= 1 .1016636 

R 2 

p = 1 .2929 

a l 

= -0.0160353 

r 3 

C T = -0.000248411 
■“■u 

a 0 

= -0.0210558 

r 4 

Sx = - 1 * 115 



r 5 

g = 9.80665 



r 6 

U ss = 77.12 



r 7 

Yss = "0.052359878 



r 8 

O u = 2.0 



r 9 

cr w = 0.0 



s 0 

C % = °' 529 

a 21 

= 5.9757330 

S] 

C kx = 4 * 87 

a 23 

= 55.0128915 

S2 

C L £ = 0.283 

a 32 

= -73.9231523 

s 3 

C L ^ = 0.0889 

a 22 

= 4.2010871 

S 4 

c D,o = 0.038 

a l 2 

= 0.0316972 

s 5 

0.0 

a l 3 

= 0.2546699 

s 6 

0.0 

a 26 

= 0.8064467 

s? 

- -0-241 

a 25 

= 0.6856605 

S8 

CmQ = -0*707 

a 35 

= 0.5113523 

s 9 

^m,o = 0.0 

a l 4 

= 0.0007159 

R A 

c = 7.0104 

a 31 

= -9.7126460 

r B 

S = 267.1 

a 33 

= 1 .5369077 


A = 7.03 

a 4 : 

= 78.1242394 

r d 

C Lf0 = 0.705 


r e 

0.0 



I 

20 
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Check Case for Program 2 


Register 

Input storage 

Output 

R 0 

g = 9.80665 

0.0 

R 1 

U ss = 77.12 

0.6027688 

r 2 

p = 1 .2929 

0. 30891 44 

r 3 

S = 267.1 

0.0064130 

r 4 

b = 43.4 

-11.394069 

r 5 

k x = 6.559296 


r 6 

a u = 2.0 

bi 5 - -0.1 779356 

r 7 

a w = -0.5 

b 1 7 = 0.0473607 

r 8 

k xz = -1 .2801 6 

b 42 = 0.0380903 

r 9 

k z = 1 2.24991 2 

b 43 = 0.01 0921 0 

s 0 

C n = -0.057 
n r 

b 24 = 0.1862234 

Si 

C n6 = 0.173 

b 35 = -0.5652042 

s 2 

C n = -0.0182 
n P 

bi 6 = 0.0594608 

S 3 

= 

b 34 = 0.0 

S 4 

C 2g = -0.21 

b 33 = 2.3929305 

s 5 

c 7 = -o.m 
*p 

b-| 4 = 1 .2648347 

s 6 

C 7 = 0.0614 

b 23 = -0.6996473 

s 7 

C Z3 = 0,0 

bs2 = 0.0 

s 8 

Cy^ = —0.866 

b 31 = 0.1 268488 

s 9 

Cy^ = 0.0 

C 30 =1.0 

r a 

m = 90909.1 

B 22 = -0.0001293 

r b 

Cy r = 0.0881 

C 21 = 0.9870954 


Y ss = -0.052359878 

b 13 = -0.1278111 

R D 

C Y = 0.0539 
y P 

Cn = -0.0078951 

r e 

0.0 


I 

24 

21 
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Check Case for Program 3 - Label A 


Register 

Input storage 


R 0 

O 

• 

o 

a 4 

R 1 

0.68541 41 

a 3 

R 2 

0.3031 611 

a 2 

r 3 

0.006391 5 

a l 

r 4 

-490.2586228 

a 0 

r 5 



r 6 

b 1 5 = -0.1779356 


*7 

bi 7 = 0.0473607 


R 8 

b 42 = 0.0380903 


r 9 

b 43 = 0.01 0921 0 


so 

b 24 = 0.1862234 


Si 

b 35 = -0.5652042 


s 2 

bi <r = 0.0594608 

1 VI 


S3 

b 34 = 0.0 


S 4 

b 33 = 2.3929305 


S5 

b! 4 = 1.2648347 


S6 

b 23 = -0.6996473 


S7 

b32 = 0.0 


S8 

b 31 = 0.1 268488 


S 9 

C 30 = 1.0 


r A 

b 2 2 = -0.0110087 


r B 

C 21 = 0.9870954 


Rc 

b 1 3 = -0.1 27381 4 


r d 

Cn = -0.0421 244 


Rg 



I 

21 



a They are normalized coefficients for the quintic; 


Output 3 

= 1 .5838890 
= 0.9679675 
= 1 .1621140 
= 0.0095553 
= -0.0001405 


thus, as = 1.0. 
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Check Case for Program 3 - Label B 


Output 


Register 

Input 

(coefficient 
of quartic) 

Ro 

a 4 = 1 .583889 

a 3 = 1 .5915011 

R 1 

a 3 = 0.9679675 

a 2 = 0.9800821 

R2 

a 2 = 1.1621140 

a] = 1 .1695744 

r 3 

a 1 = 0.0095552 

a 0 = 0.0184581 

r 4 

a 0 = -0.0001405 


*5 

Re 

R7 

Rs 

r 9 

s o 

Sl 

s 2 

s 3 

S 4 

s 5 

s 6 

s? 

s 8 

s 9 

Ra 

Rb 

Rc 

r d 

Re 


The stack contains the quintic data as follows: 
Stack register T Number of iterations 

Stack register Z tp or t-| / 2 


12 

tp = 91.0398496 (displayed as 
negative number) 


Stack register Y Root 

Stack register X 8.0 


0.0076121 

Indicates root has been found 


Use the Rl- to move data into the X register for recording. 
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Check Case for Program 4 


Store: 


a 3 = 1 

.4007102 in Rq 


a 2 ~ 1 

.1058038 in Rt 


a l 35 - 

0.0158317 in R 2 


a 0 55 ~ 

0.0227494 in R 3 


Results: 

R 0 

Root indicator 

-1.00 (indicates complex roots) 

R 1 

Real part 

-0.6946683 

r 2 

Imaginary part 

0.7924165 

s 0 

Root indicator 

0.00 (indicates real roots) 

S] 

First real root 

-0.1489289 

s 2 

Second real root 

0.1375553 


53 



APPENDIX G 


Check Case 


Register 

Input 

Ro 

-1 .0 

R 1 

-0.6946683 

r 2 

0.7924165 

r 3 


r 4 


r 5 


Re 


r 7 


Rs 


r 9 



s 0 

0 

• 

0 

S 1 

-0.1489288 

s 2 

0.1375553 

S 3 


s 4 



s 5 

S 6 

s 7 

S 8 

s 9 

Ra 

r b 

Rc 

Rd 

r e 

I 


for Program 5 

Output 


-1.0 Root 

-0 . 6946683> indicator 
0.7924165J and roots 

t-,/2 = 0.9975984 
A = 5.5228662 
N-) /2 = 0.1254783 
t = 7.9291450 
(0 n = 1 .0537969 
£ = 0.6592051 


0.0 > Root 

-0.1489288> indicator 
0.1375553J and roots 

ti /2 = 4.6532303 First root 
t D = -5.0379738 Second root 
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Check Case for Program 6 


Space axes (X S p, Y S p, Z S p) to body axes 
x sp = Ysp “ z sp = ** • ® 

Results: 

Xfo - 1 . 1 351 in R 5 
Y 5 = 1 . 0267 in R 7 
Zfc = 0.8109 in Rq 


( x b' Y b' z b) • 

y\) = 25°; 6 = 10°; $ = 30° 


Body axes (Xb/Y^Zj-j) to space axes (X sp , Y sp/ Z sp ) : 

Xb = 1.1351 y b = 1.0267 z b = 0.8109 = 25°; 0 = 10°; <t> = 30° 

Results: 

x sp = 1 • 0000 i n r 6 

*» — i nonn i n Dn 

sp 1 ' wvvyv/ -* / 

z sp = ^ • 0000 in Rg 
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A DISCUSSION OF FERRARI' S METHOD FOR THE SOLUTION OF A QUARTIC EQUATION 


Ferrari (1522-1575), an Italian mathematician, obtained the solution of a 
quartic by reducing the problem to the solution of two quadratic equations. As 
the details of obtaining the quadratic equations are not consistent among 
authors, the details of obtaining the quadratics used for the solution in this 
paper are presented. 


The general quartic equation is 

x^ + a 3 X^ + a 2 X^ + a-jx + ao = 0 
Rewrite this equation as 

x^ + a3X^ = -a2X^ - a]X - ao 
and complete the square 



(HI) 


(H2) 


(H3) 


Now, add 
variable 



to each side of equation (H3) , 


y being a dummy 



The left-hand side of equation (H4) is a perfect square. If the right-hand 
side is also a perfect square, it can be written as the square of a linear func- 
tion of x, say Cx + D. Thus, the pair of quadratics that must be solved for 
the roots of the quartic are 


x2 


a 3 y 

+ — x + — 
2 2 


± (Cx + D) 


(H5) 


The right-hand side of equation (H4) is a perfect square if, 
its discriminant is 0 


and only if. 
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= 0 


(H6) 


In this equation y has not been defined, and if equation (H6) is written as 
a function of y, it becomes 


,3 - 


a 2 y 2 + (a 3 a*i - 4a 0 )y + 


a 0 (4a 2 - 


a 3 2 ) - 


a V 


= 0 


(H7) 


This equation is called the resolvent cubic and any root y^ of equation (H7) 
insures that equation (H6) is 0. 

All that remains is the determination of the coefficients C and D. The 
discriminant equation (H6) 



permits the right-hand side of equation (H4) 


to be written as 



which is a perfect square, and the coefficients C and D are 
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i 


only if D ^ 0. The right-hand side of equation (H4) as written is a 
perfect square because 



C 



y 


and are used in place of equation (H8) if D = 0. 


/ 


(H10) 
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Stability axes, X, Y, Z 
and principal body axes 

S' S' S 


sp 

Space axes, X , Y , Z 

sp sp sp 

and principal body axes, 

S' S' S 

Figure 1.- Coordinate systems and Euler angles* Order of rotation for Euler 
angles is 0, and 4>. Moving axes translate with airplane and remain 

parallel to Earth fixed axes* Positive directions are shown* 


Earth fixed axes, 
X , Y , Z 
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